# Replication files for 
# Masterson and Lehmann
# "Refugees, Mobilization, and Humanitarian Aid: Evidence from the Syrian Refugee Crisis in Lebanon"

# File 1: Main Paper Analysis

####
# Sec. Preample -----------------------------------------------------------

## Clear all objects
remove(list = ls())

## Set workshop directory
setwd("~/Dropbox/replication/refugees-mobilization-aid-replication-files/")

# Load data
load("data/LebanonCashData.RData")

packages <- c("rdrobust", "sandwich", "robust", "xtable", 
              "apsrtable", "data.table", "lmtest", "stargazer", 
              "rdd", "texreg", "FactoMineR", "randomizr", "aod", 
              "scales", "multiwayvcov")
lapply(packages, require, character.only = T)


# Define cluster-robust standard error function  
clx <-
  function(fm, dfcw, cluster){
    library(sandwich)
    library(lmtest)
    M <- length(unique(cluster))
    N <- length(cluster)
    dfc <- (M/(M-1))*((N-1)/(N-fm$rank))
    u  <- apply(estfun(fm),2,
                function(x) tapply(x, cluster, sum))
    vcovCL <- dfc*sandwich(fm, meat=crossprod(u)/N)*(fm$df / (fm$df - (M -1)))
    coeftest(fm, vcovCL) }


# Define Cohen's d function

cohensD <- function(x, z){
  x1_bar <- mean(x[z==1])
  x2_bar <- mean(x[z==0])
  print((x1_bar-x2_bar)/sd(x))
}



# Define Stars wars function
star.wars <- function(x){
  out <- ifelse(x <= 0.1, ifelse(x <= 0.05, ifelse(x <= 0.01, "***", "**"), '*'), "")
  out
}


# Define standardization function
standardize <- function(variable){
  demeaned <- variable - mean(na.omit(variable))
  sd <- sqrt(var(na.omit(variable)))
  return(demeaned/sd)
}


##########
# Definitions -------------------------------------------------------------


# Define long outcomes data -----------------------------------------------

numVar_short <- 5
numVar_long <- 10
numVar_short2 <- numVar_long-numVar_short
numVar_PCA <- 3
numVar_total <- numVar_long + numVar_PCA

Treatment <- outcomesData_short[, "Treatment"]

write(numVar_total, file = "output/numVar_total.txt")



##########
# Define holder vectors for coefficients and standard errors
coefs_linear <- c()
coefs_linear_controls <- c()
coefs_quadratic <- c()
coefs_quadratic_controls <- c()


sesCR_linear <- c()
sesCR_linear_controls <- c()
sesCR_quadratic <- c()
sesCR_quadratic_controls <- c()


##########
# Define holder vectors for balance test 
# coefficients and standard errors
balance_coefs_linear <- c()
balance_coefs_linear_controls <- c()
balance_coefs_quadratic <- c()
balance_coefs_quadratic_controls <- c()


balance_sesCR_linear <- c()
balance_sesCR_linear_controls <- c()
balance_sesCR_quadratic <- c()
balance_sesCR_quadratic_controls <- c()






##########
# Sharing tests: Define holder vectors for coefficients and standard errors
coefs_linear53_37 <- c()
coefs_linear_controls53_37 <- c()
coefs_quadratic53_37 <- c()
coefs_quadratic_controls53_37 <- c()


sesCR_linear53_37 <- c()
sesCR_linear_controls53_37 <- c()
sesCR_quadratic53_37 <- c()
sesCR_quadratic_controls53_37 <- c()


##########
# Sharing tests: Define holder vectors for balance test coefficients and standard errors
balance_coefs_linear53_37 <- c()
balance_coefs_linear_controls53_37 <- c()
balance_coefs_quadratic53_37 <- c()
balance_coefs_quadratic_controls53_37 <- c()


balance_sesCR_linear53_37 <- c()
balance_sesCR_linear_controls53_37 <- c()
balance_sesCR_quadratic53_37 <- c()
balance_sesCR_quadratic_controls53_37 <- c()








##########
#### Define descriptive statistics
### Output basic descriptive statistics for households, across altitudes
desRows <- c("At least one person returned to Syria (yes=1, no=0)",
             "Someone moved to Syria for money (yes=1, no=0)",
             "Number of returnees",
             "Among HHs with a return, change in number of men aged 18 to 50",
             "Number of men aged 18 to 50",
             "Dangerous work (yes=1, no=0)",
            # "Return to Syria for money (yes=1, no=0)",
            # "Number of returns",
             "Under siege (yes=1, no=0)",
             "Remittances from Syria (yes=1, no=0)"
)

desCols <- c("Mean", "Median", "SD", "Range")

mobilization_descriptives <- matrix(NA, nrow = length(desRows), ncol = length(desCols))

rownames(mobilization_descriptives) <- desRows

colnames(mobilization_descriptives) <- desCols

mobDesData <- data.frame(cbind(  
  # The share of households that had someone go back
  data[which(data$Treatment==0), Q64_0_1], 
  # Since November, have any members of the household gone back to Syria to earn money?
  data[which(data$Treatment==0), Q91_1],  
  # Average number of returnees per HH?
  data[which(data$Treatment==0), Q64],   
  # Among families that had someone go back in the last four months, 
  # What is the mean change in the nuber of men?
  (I(data[which(data$Treatment==0), rescaled_men_18_50] * data[which(data$Treatment==0), Q64_0_1])),
  data[which(data$Treatment==0), SumFightingAgeMale18_50],   # men 18-50"
  data[which(data$Treatment==0), Q89], #Dangerous work"
  data[which(data$Treatment==0), Q66b1],     #   Under siege"
  as.numeric(data[which(data$Treatment==0), Q51_1]>0)  #"Remittances from Syria"
))

colnames(mobDesData) <- desRows



#####   Input content
for(i in 1:dim(mobDesData)[2]){
  
  mobilization_descriptives[i, "Mean"] <- round(mean(mobDesData[, i]), digits = 3)
  
  mobilization_descriptives[i, "Median"] <- round(median(mobDesData[, i]), digits = 3)
  
  mobilization_descriptives[i, "SD"] <- round(sd(mobDesData[, i]), digits = 2)
  
  minLoopHolder <- round(min(mobDesData[, i]), digits = 3)
  maxLoopHolder <- round(max(mobDesData[, i]), digits = 3)
  
  mobilization_descriptives[i, "Range"] <- paste(minLoopHolder, maxLoopHolder, sep = "--")
  
}



#####
## Output results
holder <- colnames(mobilization_descriptives)
mobilization_descriptives <- cbind(rownames(mobilization_descriptives), mobilization_descriptives)
rownames(mobilization_descriptives) <- c()
colnames(mobilization_descriptives) <- c("Outcome", holder)

print(xtable(data.frame(mobilization_descriptives), label = "mobDescriptives"), 
      only.contents=TRUE, include.rownames=F, 
      include.colnames=T, floating=F,
      #hline.after=numVar_short*3, 
      #	size="\\fontsize{9pt}{10pt}\\selectfont",
      file = "output/mobilizationDescriptivesTable.tex")






#############################################
# Results ---------------------------------------------------
####   Balance Results  
####        Balance Tests 
## Define balance variables matrix
## Run balance tests



# Define a treatment indicator vector
Treatment_balance <- balance_variables[, "Treatment"]

# Define a variable for the number of balance tests, subtracting the three non-balance test variables included in the matrix
numBalanceTests <- dim(balance_variables)[2]-3



#############################
####  Linear balance analysis -----------------------------------------------
##### 	 balance without controls
####### 	   balance clx linear no controls

for(i in 1:numBalanceTests){
  # Drop rows with an NA for the variable being balance tested
  balance_variables_NA_drop <- balance_variables[which(is.na(balance_variables[, i])==0), ]
  #Define the lm model
  modelLinNoCont <- lm(balance_variables_NA_drop[, i]  ~ balance_variables_NA_drop[, "Treatment"] + balance_variables_NA_drop[, "WINTELEV_oct13_rescaled"] +  I(balance_variables_NA_drop[, "WINTELEV_oct13_rescaled"] * balance_variables_NA_drop[, "Treatment"])  )
  # 	Optional: print the colnames to follow the regressions
  # 	print(colnames(balance_variables)[i])
  #	Optional: print the regression results to follow  
  #	print(clx(fm = modelLinNoCont,  cluster = balance_variables_NA_drop$cas_code))
  # Add the calculated coefficient to the coefs results matrix
  balance_coefs_linear[i] <- clx(fm = modelLinNoCont,  cluster = balance_variables_NA_drop$cas_code)[2, 1]
  # Add the calculated cluster-robust standard error to the CR SEs results matrix
  balance_sesCR_linear[i] <- clx(fm = modelLinNoCont,  cluster = balance_variables_NA_drop$cas_code)[2, 2]
}



####################
####   quadratic balance analysis ------------------------------------------
#####    balance without controls 
####### 	  balance clx quadratic no controls

for(i in 1:numBalanceTests){
  # Drop rows with an NA for the variable being balance tested
  balance_variables_NA_drop <- balance_variables[which(is.na(balance_variables[, i])==0), ]
  modelQuadNoCont <- lm(balance_variables_NA_drop[, i]  ~ balance_variables_NA_drop[, "Treatment"] + balance_variables_NA_drop[, "WINTELEV_oct13_rescaled"] +  I(balance_variables_NA_drop[, "WINTELEV_oct13_rescaled"] * balance_variables_NA_drop[, "Treatment"]) +  balance_variables_NA_drop[, "WINTELEV_oct13_rescaled"]^2  +  I(balance_variables_NA_drop[, "Treatment"] * balance_variables_NA_drop[, "WINTELEV_oct13_rescaled"]^2) )
  #	print(colnames(balance_variables)[i])
  #	print(clx(fm = modelQuadNoCont,  cluster = balance_variables_NA_drop$cas_code))
  balance_coefs_quadratic[i] <- clx(fm = modelQuadNoCont,  cluster = balance_variables_NA_drop$cas_code)[2, 1]
  balance_sesCR_quadratic[i] <- clx(fm = modelQuadNoCont,  cluster = balance_variables_NA_drop$cas_code)[2, 2]
}



p_valuesCR_balance <- matrix(NA, nrow = numBalanceTests, ncol = 2)
colnames(p_valuesCR_balance) <- c("Linear", "Quadratic")



#bind the results into matrices
coefs_balance <- matrix(cbind(balance_coefs_linear, balance_coefs_quadratic), ncol = 2)
colnames(coefs_balance) <- c("Linear", "Quadratic")

sesCR_balance <- matrix(cbind(balance_sesCR_linear, balance_sesCR_quadratic), ncol = 2)
colnames(sesCR_balance) <- c("Linear", "Quadratic")



# generate t-stat and p values
for(i in 1:2){
  
  zCR_balance <- coefs_balance[, i] / sesCR_balance[, i]  # Calculate Z-score for each coef
  p_valuesCR_balance[, i] <- 2*pnorm(-abs(zCR_balance))  # Calculate p-value for each coef
 
}


# Review output
sum(p_valuesCR_balance[, "Linear"] <=0.01, na.rm=T)/numBalanceTests
sum(p_valuesCR_balance[, "Quadratic"] <=0.01, na.rm=T)/numBalanceTests

sum(p_valuesCR_balance[, "Linear"] <=0.05, na.rm=T)/numBalanceTests
  write(round(sum(p_valuesCR_balance[, "Linear"] <=0.05, na.rm=T)/numBalanceTests, 2)*100, file = "output/percent_imbalance_5percent_linear.txt")
sum(p_valuesCR_balance[, "Quadratic"] <=0.05, na.rm=T)/numBalanceTests
  write(round(sum(p_valuesCR_balance[, "Quadratic"] <=0.05, na.rm=T)/numBalanceTests, 2)*100, file = "output/percent_imbalance_5percent_quadratic.txt")

sum(p_valuesCR_balance[, "Linear"] <=0.10, na.rm=T)/numBalanceTests
sum(p_valuesCR_balance[, "Quadratic"] <=0.10, na.rm=T)/numBalanceTests



# Review output
which(p_valuesCR_balance[, "Linear"] <=0.05)
which(p_valuesCR_balance[, "Quadratic"] <=0.05)
which(p_valuesCR_balance[, "Linear"] <=0.05 & p_valuesCR_balance[, "Quadratic"] <=0.05)


write(round(length(which(p_valuesCR_balance[, "Linear"] <=0.05 & p_valuesCR_balance[, "Quadratic"] <=0.05)), 2), file = "output/number_imbalance_5percent_both_linear_and_quadratic.txt")
write(round(length(which(p_valuesCR_balance[, "Linear"] <=0.05 & p_valuesCR_balance[, "Quadratic"] <=0.05))/numBalanceTests, 2)*100, file = "output/percent_imbalance_5percent_both_linear_and_quadratic.txt")


colnames(balance_variables)[which(p_valuesCR_balance[, "Linear"] <=0.05)]
colnames(balance_variables)[which(p_valuesCR_balance[, "Quadratic"] <=0.05)]
colnames(balance_variables)[which(p_valuesCR_balance[, "Linear"] <=0.05 & p_valuesCR_balance[, "Quadratic"] <=0.05)]


sig_imbalance <- 
  cbind(c( 
    (sum(p_valuesCR_balance[, "Linear"] <=0.01, na.rm=T)/numBalanceTests), 
    (sum(p_valuesCR_balance[, "Linear"] <=0.05, na.rm=T)/numBalanceTests),  
    (sum(p_valuesCR_balance[, "Linear"] <=0.10, na.rm=T)/numBalanceTests)  ),
  c(  (sum(p_valuesCR_balance[, "Quadratic"] <=0.01, na.rm=T)/numBalanceTests),
      (sum(p_valuesCR_balance[, "Quadratic"] <=0.05, na.rm=T)/numBalanceTests),
      (sum(p_valuesCR_balance[, "Quadratic"] <=0.10, na.rm=T)/numBalanceTests)  ))

colnames(sig_imbalance) <- c("Linear", "Quadratic")
rownames(sig_imbalance) <- c("p<0.01", "p<0.05", "p<0.10")

print(xtable(sig_imbalance), only.contents=TRUE, include.rownames=T, include.colnames=T, floating=F, file = "output/balance_imbalance.tex")


# which variables are imbalanced in both specifications
colnames(balance_variables)[which(p_valuesCR_balance[, "Linear"] <=0.05)[which(p_valuesCR_balance[, "Linear"] <=0.05) %in% which(p_valuesCR_balance[, "Quadratic"] <=0.05)]]

write(numBalanceTests,  file = "output/numBalanceTests.txt")
write(sum(p_valuesCR_balance[, "Linear"] <=0.05, na.rm=T),  file = "output/numImbalanced.txt")

# The sign of the imbalance for the imbalanced vars
balance_coefs_quadratic[which(p_valuesCR_balance[, "Linear"] <=0.05)[which(p_valuesCR_balance[, "Linear"] <=0.05) %in% which(p_valuesCR_balance[, "Quadratic"] <=0.05)]]

# imbalance in year
year_imbalance <- abs(c(round(balance_coefs_linear[which(colnames(balance_variables)=="Q11_Year")], digits = 2),
round(balance_coefs_quadratic[which(colnames(balance_variables)=="Q11_Year")], digits = 2)))

write(min(year_imbalance), file = "output/min_year_imbalance.txt")
write(max(year_imbalance), file = "output/max_year_imbalance.txt")

balance_coefs_linear[which(colnames(balance_variables)=="SSFSC")]
balance_coefs_quadratic[which(colnames(balance_variables)=="SSFSC")]



####
# Numbers for Main Paper -------------------------------------------


gendersHolder <- (na.omit(c(data$Q64a_1_1, data$Q64a_1_2, 
                            data$Q64a_1_3, data$Q64a_1_4, 
                            data$Q64a_1_5, data$Q64a_1_6,
                            data$Q64a_1_7, data$Q64a_1_8,
                            data$Q64a_1_9, data$Q64a_1_10   )))


gendersHolder_control <- (na.omit(c(data[which(data$Treatment==0), Q64a_1_1], 
                                    data[which(data$Treatment==0), Q64a_1_2], 
                                    data[which(data$Treatment==0), Q64a_1_3], 
                                    data[which(data$Treatment==0), Q64a_1_4], 
                                    data[which(data$Treatment==0), Q64a_1_5], 
                                    data[which(data$Treatment==0), Q64a_1_6],
                                    data[which(data$Treatment==0), Q64a_1_7], 
                                    data[which(data$Treatment==0), Q64a_1_8],
                                    data[which(data$Treatment==0), Q64a_1_9], 
                                    data[which(data$Treatment==0), Q64a_1_10]   )))



#################################
###   Save average household labor income
write(round(mean(data$Q18_USD, na.rm=TRUE), digits = 2), file = "output/Q18_USD_mean_labor_income.txt")
write(round(mean(data$Q18_USD[Treatment==0], na.rm=TRUE)), file = "output/Q18_USD_mean_labor_income_control.txt")

#   controlGroupSize.txt
# How many control group observations are there?
write(sum(data$Treatment==0), file = "output/controlGroupSize.txt")

#  TotalHouseholdsWithMoversToSyria_control.txt
write( sum(is.na( data[which(data$Treatment==0), Q64a_2_1])==0), 
       file = "output/TotalHouseholdsWithMoversToSyria_control.txt")

# How many movers to Syria were there total?
write(length(gendersHolder), 
      file = "output/TotalMoversToSyria.txt")

#  ShareHouseholdsWithMoversToSyria_control.txt
write(length(gendersHolder_control), 
      file = "output/TotalMoversToSyria_control.txt")

#  AvgNumReturnsAmongReturns_control.txt
write( round(   mean(data[which(data$Q64>0 & data$Treatment==0), Q64])  , 1), 
       file = "output/AvgNumReturnsAmongReturns_control.txt")

#  ShareHouseholdsWithMoversToSyria_control_decimal.txt
write( round((sum(is.na( data[which(data$Treatment==0), Q64a_2_1])==0) / dim( data[which(data$Treatment==0), ])[1]), 3), 
       file = "output/ShareHouseholdsWithMoversToSyria_control_decimal.txt")

#  AvgNumReturnsAmongReturns_control.txt
write( round(   mean(data[which(data$Q64>0 & data$Treatment==0), Q64])  , 1), 
       file = "output/AvgNumReturnsAmongReturns_control.txt")

# ShareHouseholdsWithMoversToSyria_control.txt
write( round((sum(is.na( data[which(data$Treatment==0), Q64a_2_1])==0) / dim( data[which(data$Treatment==0), ])[1]), 3)*100, 
       file = "output/ShareHouseholdsWithMoversToSyria_control.txt")

#  AvgNumReturnsAmongReturns_control.txt
write( round(   mean(data[which(data$Q64>0 & data$Treatment==0), Q64])  , 1), 
       file = "output/AvgNumReturnsAmongReturns_control.txt")




#
#
#
#
#
#


##############################
#### Define descriptive statistics ---------------------------------------
### Output basic descriptive statistics for households, across altitudes
desStatsCols <- c("Mean", "Median", "SD", "Range")

desStatsRows_short <- c( "DiD18_50",
                         "Q64>0",
                         "I(data$data$rescaled_men_18_50 * data$Q64_0_1)",
                         "66b1",
                         "returnedXundersiege") 

# Define Variable Means for Control Group ----------------------------------
basicStats_short_control <- matrix(NA, nrow = length(desStatsRows_short), ncol = length(desStatsCols))
rownames(basicStats_short_control) <- desStatsRows_short
colnames(basicStats_short_control) <- desStatsCols

#####   Input content -- descriptive stats for key short outcomes
for(i in 1:numVar_short){
  basicStats_short_control[i, "Mean"] <- round(mean(outcomesData_short[which(data$Treatment==0), i]), digits = 3)
  basicStats_short_control[i, "Median"] <- round(median(outcomesData_short[which(data$Treatment==0), i]), digits = 3)
  basicStats_short_control[i, "SD"] <- round(sd(outcomesData_short[which(data$Treatment==0), i]), digits = 2)
  minLoopHolder <- round(min(outcomesData_short[which(data$Treatment==0), i]), digits = 3)
  maxLoopHolder <- round(max(outcomesData_short[which(data$Treatment==0), i]), digits = 3)
  basicStats_short_control[i, "Range"] <- paste(minLoopHolder, maxLoopHolder, sep = "--")
}



####
# Main RD regressions -----------------------------------------------------

####
#   All ten outcomes and 3 PCA dimensions


###############################
# Define holder vectors for coefficients and standard errors
coefs_long_linear <- c()
coefs_long_linear_controls <- c()
coefs_long_quadratic <- c()
coefs_long_quadratic_controls <- c()


sesCR_long_linear <- c()
sesCR_long_linear_controls <- c()
sesCR_long_quadratic <- c()
sesCR_long_quadratic_controls <- c()





##############################
# results_long --------------------------------------------

####   Linear analysis
##### 	   without controls
####### 	  clx linear no controls

for(i in 1:numVar_total){
  modelLinNoCont_long <- lm(outcomesData_long[, i]  ~ outcomesData_long[, "Treatment"] + outcomesData_long[, "WINTELEV_oct13_rescaled"] +  I(outcomesData_long[, "WINTELEV_oct13_rescaled"] * outcomesData_long[, "Treatment"])  )

  coefs_long_linear[i] <- clx(fm = modelLinNoCont_long,  cluster = outcomesData_long$cas_code)[2, 1]
  sesCR_long_linear[i] <- clx(fm = modelLinNoCont_long,  cluster = outcomesData_long$cas_code)[2, 2]
}


##################################################
##### 	with controls
#######    clx linear with controls

for(i in 1:numVar_total){
  modelLinWCont_long <- lm(outcomesData_long[, i]  ~ outcomesData_long[, "Treatment"] + outcomesData_long[, "WINTELEV_oct13_rescaled"] +  I(outcomesData_long[, "WINTELEV_oct13_rescaled"] * outcomesData_long[, "Treatment"]) +
                             outcomesData_long$age0to4 +
                             outcomesData_long$age5to12 +
                             outcomesData_long$age13to17 +
                             outcomesData_long$age18to59 +
                             outcomesData_long$age60plus +
                             outcomesData_long$Oct2013Pop +
                             outcomesData_long$disabled +
                             outcomesData_long$origin_gov_0 +
                             outcomesData_long$origin_gov_4 +
                             outcomesData_long$origin_gov_6 +
                             outcomesData_long$origin_gov_10 +
                             outcomesData_long$educ_le_priinc +
                             outcomesData_long$educ_pricom +
                             outcomesData_long$educ_middle +
                             outcomesData_long$educ_sec_tech_uni +
                             outcomesData_long$age_head +
                             outcomesData_long$monthinleb +
                             outcomesData_long$base_relatives +
                             outcomesData_long$base_friends 
  )

  
  coefs_long_linear_controls[i] <- clx(fm = modelLinWCont_long,  cluster = outcomesData_long$cas_code)[2, 1]
  sesCR_long_linear_controls[i] <- clx(fm = modelLinWCont_long,  cluster = outcomesData_long$cas_code)[2, 2]

    
}



############################################################
####  quadratic analysis
##### 	 without controls 
####### 	  clx quadratic no controls

for(i in 1:numVar_total){
  modelQuadNoCont_long <- lm(outcomesData_long[, i]  ~ outcomesData_long[, "Treatment"] + outcomesData_long[, "WINTELEV_oct13_rescaled"] +  I(outcomesData_long[, "WINTELEV_oct13_rescaled"] * outcomesData_long[, "Treatment"]) +  outcomesData_long[, "WINTELEV_oct13_rescaled"]^2  +  I(outcomesData_long[, "Treatment"] * outcomesData_long[, "WINTELEV_oct13_rescaled"]^2) )
  #print(colnames(outcomesData_long)[i])
  #print(clx(fm = modelQuadNoCont_long,  cluster = outcomesData_long$cas_code))
  coefs_long_quadratic[i] <- clx(fm = modelQuadNoCont_long,  cluster = outcomesData_long$cas_code)[2, 1]
  sesCR_long_quadratic[i] <- clx(fm = modelQuadNoCont_long,  cluster = outcomesData_long$cas_code)[2, 2]
}







###############################
##### 	 ith controls
####### 	   clx quadratic with controls
for(i in 1:numVar_total){
  modelQuadWCont_long <- lm(outcomesData_long[, i]  ~ outcomesData_long[, "Treatment"] + outcomesData_long[, "WINTELEV_oct13_rescaled"] +  I(outcomesData_long[, "WINTELEV_oct13_rescaled"] * outcomesData_long[, "Treatment"]) + outcomesData_long[, "WINTELEV_oct13_rescaled"]^2  +  I(outcomesData_long[, "Treatment"] * outcomesData_long[, "WINTELEV_oct13_rescaled"]^2) +
                              outcomesData_long$age0to4 +
                              outcomesData_long$age5to12 +
                              outcomesData_long$age13to17 +
                              outcomesData_long$age18to59 +
                              outcomesData_long$age60plus +
                              outcomesData_long$Oct2013Pop +
                              outcomesData_long$disabled +
                              outcomesData_long$origin_gov_0 +
                              outcomesData_long$origin_gov_4 +
                              outcomesData_long$origin_gov_6 +
                              outcomesData_long$origin_gov_10 +
                              outcomesData_long$educ_le_priinc +
                              outcomesData_long$educ_pricom +
                              outcomesData_long$educ_middle +
                              outcomesData_long$educ_sec_tech_uni +
                              outcomesData_long$age_head +
                              outcomesData_long$monthinleb +
                              outcomesData_long$base_relatives +
                              outcomesData_long$base_friends 
  )
  #print(colnames(outcomesData_long)[i])
  #print(clx(fm = modelQuadWCont_long,  cluster = outcomesData_long$cas_code))
  coefs_long_quadratic_controls[i] <- clx(fm = modelQuadWCont_long,  cluster = outcomesData_long$cas_code)[2, 1]
  sesCR_long_quadratic_controls[i] <- clx(fm = modelQuadWCont_long,  cluster = outcomesData_long$cas_code)[2, 2]
}







#bind the results_long into matrices

coefs_long_unrounded <- matrix(cbind(coefs_long_linear, coefs_long_linear_controls, coefs_long_quadratic, coefs_long_quadratic_controls), ncol = 4)
sesCR_long_unrounded <- matrix(cbind(sesCR_long_linear, sesCR_long_linear_controls, sesCR_long_quadratic, sesCR_long_quadratic_controls), ncol = 4)







#
#
#
#
#
#
#
#
#
#
#
####  Output results_long table for primary results_long


# Four models, 10 variables each. 
# Define the row names of the text output
# Blank row names will be where SEs are located
# Create empty matrices to store results_long below 

# p values cluster-robust
p_valuesCR_long <- matrix(NA, nrow = numVar_total, ncol = 4)



# Define Names
colnames(coefs_long_unrounded) <- c("Linear", "Linear_C", "Quadratic", "Quadratic_C")
rownames(coefs_long_unrounded) <- as.character(paste("V", 1:numVar_total, sep = ""))

colnames(sesCR_long_unrounded) <- c("Linear", "Linear_C", "Quadratic", "Quadratic_C")
rownames(sesCR_long_unrounded) <- as.character(paste("V", 1:numVar_total, sep = ""))

colnames(p_valuesCR_long) <- c("Linear", "Linear_C", "Quadratic", "Quadratic_C")
rownames(p_valuesCR_long) <- as.character(paste("V", 1:numVar_total, sep = ""))






# Definte coefficient matrix 
coefs_long_unrounded[, "Linear"] <- coefs_long_linear
coefs_long_unrounded[, "Linear_C"] <- coefs_long_linear_controls
coefs_long_unrounded[, "Quadratic"] <- coefs_long_quadratic
coefs_long_unrounded[, "Quadratic_C"] <- coefs_long_quadratic_controls

# Definte cluster-robust SE matrix 
sesCR_long_unrounded[, "Linear"] <- sesCR_long_linear
sesCR_long_unrounded[, "Linear_C"] <- sesCR_long_linear_controls
sesCR_long_unrounded[, "Quadratic"] <- sesCR_long_quadratic
sesCR_long_unrounded[, "Quadratic_C"] <- sesCR_long_quadratic_controls


# generate t-stat and p values
for(i in 1:4){
  
  zCR <- coefs_long_unrounded[, i] / sesCR_long_unrounded[, i]  # Calculate Z-score for each coef
  p_valuesCR_long[, i] <- 2*pnorm(-abs(zCR))  # Calculate p-value for each coef
  
}

## Define rounded coefficients vector including scietific notation for the one tiny negative result
coefs_long_w_text <- 
  ifelse(
    round(coefs_long_unrounded, digits = 3) == 0,
    yes = formatC(coefs_long_unrounded, format = "e", digits = 0),
    no = round(coefs_long_unrounded, digits = 3)
  )

### Round variables for display in table post significance computation
# These rounded values are not used in significance calculations
# Only the p-values calculated above are used to calculate significance
coefs_long <- round(coefs_long_unrounded, digits = 3)  
sesCR_long <- round(sesCR_long_unrounded, digits = 3)



#### Define long summary stats objects
desStatsRows_long <- c( "DiD18_50",
                        "Q64>0",
                        "I(rescaled_men_18_50xQ64_0_1)",
                        "66b1",
                        "returnedXundersiege",
                        "SumFightingAgeMale18_50",
                        "Q89",
                        "Q91_1",                    
                        "Q51_1>0",      
                        "I(rescaled_men_18_50xQ66b1)") 

#    DEFINE the LONG basics stats holder
basicStats_long <- matrix(NA, nrow = length(desStatsRows_long), ncol = length(desStatsCols))

rownames(basicStats_long) <- desStatsRows_long

colnames(basicStats_long) <- desStatsCols

#####   Input content
for(i in 1:numVar_long){
  
  basicStats_long[i, "Mean"] <- round(mean(outcomesData_long[, i]), digits = 3)
  
  basicStats_long[i, "Median"] <- round(median(outcomesData_long[, i]), digits = 3)
  
  basicStats_long[i, "SD"] <- round(sd(outcomesData_long[, i]), digits = 2)
  
  minLoopHolder <- round(min(outcomesData_long[, i]), digits = 3)
  maxLoopHolder <- round(max(outcomesData_long[, i]), digits = 3)
  
  basicStats_long[i, "Range"] <- paste(minLoopHolder, maxLoopHolder, sep = "--")
  
}

#    Define basics stats for the control group
basicStats_long_control <- matrix(NA, nrow = length(desStatsRows_long), ncol = length(desStatsCols))

rownames(basicStats_long_control) <- desStatsRows_long

colnames(basicStats_long_control) <- desStatsCols


#####   Input content -- descriptive stats for key long outcomes
for(i in 1:numVar_long){
  
  basicStats_long_control[i, "Mean"] <- round(mean(outcomesData_long[which(data$Treatment==0), i]), digits = 3)
  
  basicStats_long_control[i, "Median"] <- round(median(outcomesData_long[which(data$Treatment==0), i]), digits = 3)
  
  basicStats_long_control[i, "SD"] <- round(sd(outcomesData_long[which(data$Treatment==0), i]), digits = 2)
  
  minLoopHolder <- round(min(outcomesData_long[which(data$Treatment==0), i]), digits = 3)
  maxLoopHolder <- round(max(outcomesData_long[which(data$Treatment==0), i]), digits = 3)
  
  basicStats_long_control[i, "Range"] <- paste(minLoopHolder, maxLoopHolder, sep = "--")
  
}








####
# Regression Output for All Outcomes --------------------------------------
#  results_short1.tex 
#  results_short2.tex 


##############################
#  Define holder vectors for coefficients and standard errors
coefs_short_linear <- c()
coefs_short_linear_controls <- c()
coefs_short_quadratic <- c()
coefs_short_quadratic_controls <- c()


sesCR_short_linear <- c()
sesCR_short_linear_controls <- c()
sesCR_short_quadratic <- c()
sesCR_short_quadratic_controls <- c()





##############################
# results_short ------------------------------

####  Linear analysis
##### 	 without controls
####### 	 clx linear no controls

for(i in 1:numVar_short){
  modelLinNoCont_short <- lm(outcomesData_short[, i]  ~ outcomesData_short[, "Treatment"] + outcomesData_short[, "WINTELEV_oct13_rescaled"] +  I(outcomesData_short[, "WINTELEV_oct13_rescaled"] * outcomesData_short[, "Treatment"])  )
  #print(colnames(outcomesData_short)[i])
  #print(clx(fm = modelLinNoCont_short,  cluster = outcomesData_short$cas_code))
  coefs_short_linear[i] <- clx(fm = modelLinNoCont_short,  cluster = outcomesData_short$cas_code)[2, 1]
  sesCR_short_linear[i] <- clx(fm = modelLinNoCont_short,  cluster = outcomesData_short$cas_code)[2, 2]
}


####################
##### 	  with controls
####### 	   clx linear with controls

for(i in 1:numVar_short){
  modelLinWCont_short <- lm(outcomesData_short[, i]  ~ outcomesData_short[, "Treatment"] + outcomesData_short[, "WINTELEV_oct13_rescaled"] +  I(outcomesData_short[, "WINTELEV_oct13_rescaled"] * outcomesData_short[, "Treatment"]) +
                              outcomesData_short$age0to4 +
                              outcomesData_short$age5to12 +
                              outcomesData_short$age13to17 +
                              outcomesData_short$age18to59 +
                              outcomesData_short$age60plus +
                              outcomesData_short$Oct2013Pop +
                              outcomesData_short$disabled +
                              outcomesData_short$origin_gov_0 +
                              outcomesData_short$origin_gov_4 +
                              outcomesData_short$origin_gov_6 +
                              outcomesData_short$origin_gov_10 +
                              outcomesData_short$educ_le_priinc +
                              outcomesData_short$educ_pricom +
                              outcomesData_short$educ_middle +
                              outcomesData_short$educ_sec_tech_uni +
                              outcomesData_short$age_head +
                              outcomesData_short$monthinleb +
                              outcomesData_short$base_relatives +
                              outcomesData_short$base_friends 
  )
  #print(colnames(outcomesData_short)[i])
  #print(clx(fm = modelLinWCont_short,  cluster = outcomesData_short$cas_code))
  coefs_short_linear_controls[i] <- clx(fm = modelLinWCont_short,  cluster = outcomesData_short$cas_code)[2, 1]
  sesCR_short_linear_controls[i] <- clx(fm = modelLinWCont_short,  cluster = outcomesData_short$cas_code)[2, 2]
}



####################
####  quadratic analysis --------------------------------------------
##### 	   without controls 
#######        clx quadratic no controls

for(i in 1:numVar_short){
  modelQuadNoCont_short <- lm(outcomesData_short[, i]  ~ outcomesData_short[, "Treatment"] + outcomesData_short[, "WINTELEV_oct13_rescaled"] +  I(outcomesData_short[, "WINTELEV_oct13_rescaled"] * outcomesData_short[, "Treatment"]) +  outcomesData_short[, "WINTELEV_oct13_rescaled"]^2  +  I(outcomesData_short[, "Treatment"] * outcomesData_short[, "WINTELEV_oct13_rescaled"]^2) )
  #print(colnames(outcomesData_short)[i])
  #print(clx(fm = modelQuadNoCont_short,  cluster = outcomesData_short$cas_code))
  coefs_short_quadratic[i] <- clx(fm = modelQuadNoCont_short,  cluster = outcomesData_short$cas_code)[2, 1]
  sesCR_short_quadratic[i] <- clx(fm = modelQuadNoCont_short,  cluster = outcomesData_short$cas_code)[2, 2]
}







##############################
##### 	  with controls ---------------------------------------
####### 	     clx quadratic with controls
for(i in 1:numVar_short){
  modelQuadWCont_short <- lm(outcomesData_short[, i]  ~ outcomesData_short[, "Treatment"] + outcomesData_short[, "WINTELEV_oct13_rescaled"] +  I(outcomesData_short[, "WINTELEV_oct13_rescaled"] * outcomesData_short[, "Treatment"]) + outcomesData_short[, "WINTELEV_oct13_rescaled"]^2  +  I(outcomesData_short[, "Treatment"] * outcomesData_short[, "WINTELEV_oct13_rescaled"]^2) +
                               outcomesData_short$age0to4 +
                               outcomesData_short$age5to12 +
                               outcomesData_short$age13to17 +
                               outcomesData_short$age18to59 +
                               outcomesData_short$age60plus +
                               outcomesData_short$Oct2013Pop +
                               outcomesData_short$disabled +
                               outcomesData_short$origin_gov_0 +
                               outcomesData_short$origin_gov_4 +
                               outcomesData_short$origin_gov_6 +
                               outcomesData_short$origin_gov_10 +
                               outcomesData_short$educ_le_priinc +
                               outcomesData_short$educ_pricom +
                               outcomesData_short$educ_middle +
                               outcomesData_short$educ_sec_tech_uni +
                               outcomesData_short$age_head +
                               outcomesData_short$monthinleb +
                               outcomesData_short$base_relatives +
                               outcomesData_short$base_friends 
  )
  #print(colnames(outcomesData_short)[i])
  #print(clx(fm = modelQuadWCont_short,  cluster = outcomesData_short$cas_code))
  coefs_short_quadratic_controls[i] <- clx(fm = modelQuadWCont_short,  cluster = outcomesData_short$cas_code)[2, 1]
  sesCR_short_quadratic_controls[i] <- clx(fm = modelQuadWCont_short,  cluster = outcomesData_short$cas_code)[2, 2]
}







#bind the results_short into matrices

coefs_short <- matrix(cbind(coefs_short_linear, coefs_short_linear_controls, coefs_short_quadratic, coefs_short_quadratic_controls), ncol = 4)
sesCR_short <- matrix(cbind(sesCR_short_linear, sesCR_short_linear_controls, sesCR_short_quadratic, sesCR_short_quadratic_controls), ncol = 4)







#
#
#
#
#
#
#
#
#
#
#
####  Output results_short table for primary results_short


# Four models, twelve variables each. 

# Define a variable for the number of rows in the tex output
numRows_short <- numVar_short*3

# Define the row names of the text output
# Blank row names will be where SEs are located
RowsNameswSEs_short <- c("$\\Delta$ men 18-50", "", "", 
                         "HH member returned to Syria", "", "", 
                         "($\\Delta$ men 18-50 x returned to Syria)",  "", "", 
                         "Family member under siege (yes=1, no=0)", "", "",
                         #       "I(DiD Men 18-50 x under siege)", "", ""
                         "(Returned x under siege)", "", ""
) 



# Create empty matrices to store results_short below 


# p values cluster-robust
p_valuesCR_short <- matrix(NA, nrow = numVar_short, ncol = 4)



# Define Names
colnames(coefs_short) <- c("Linear", "Linear_C", "Quadratic", "Quadratic_C")
rownames(coefs_short) <- as.character(paste("V", 1:numVar_short, sep = ""))

colnames(sesCR_short) <- c("Linear", "Linear_C", "Quadratic", "Quadratic_C")
rownames(sesCR_short) <- as.character(paste("V", 1:numVar_short, sep = ""))

colnames(p_valuesCR_short) <- c("Linear", "Linear_C", "Quadratic", "Quadratic_C")
rownames(p_valuesCR_short) <- as.character(paste("V", 1:numVar_short, sep = ""))






# Definte coefficient matrix 
coefs_short[, "Linear"] <- coefs_short_linear
coefs_short[, "Linear_C"] <- coefs_short_linear_controls
coefs_short[, "Quadratic"] <- coefs_short_quadratic
coefs_short[, "Quadratic_C"] <- coefs_short_quadratic_controls

# Definte cluster-robust SE matrix 
sesCR_short[, "Linear"] <- sesCR_short_linear
sesCR_short[, "Linear_C"] <- sesCR_short_linear_controls
sesCR_short[, "Quadratic"] <- sesCR_short_quadratic
sesCR_short[, "Quadratic_C"] <- sesCR_short_quadratic_controls

# generate t-stat and p values
for(i in 1:4){
  
  zCR <- coefs_short[, i] / sesCR_short[, i]  # Calculate Z-score for each coef
  p_valuesCR_short[, i] <- 2*pnorm(-abs(zCR))  # Calculate p-value for each coef
  
}
# Create rounded p value matrix
p_valuesCR_short_rounded <- round(p_valuesCR_short, 2)

### Round variables for display in table post significance computation
# These rounded values are not used in significance calculations
# Only the p-values calculated above are used to calculate significance
coefs_short <- round(coefs_short, digits = 3)  
sesCR_short <- round(sesCR_short, digits = 3)


###### For appendix
# Output short1 results table for all models, first five variables
# All models
# columns are models
# rows are point estimate, standard error of the point estimate, p-value, quasi-cohen's d for beta hat, and quasi-cohen's d for beta hat confidence interval range

RowsNameswSEs_short1 <- c("$\\hat{\\beta}$: $\\Delta$ men 18-50", "", "  \\emph{p}-value", "  Standarized $\\hat{\\beta}$", "Standardized 95\\% CI range", 
                          "$\\hat{\\beta}$: HH member returned to Syria", "", "  \\emph{p}-value", "  Standarized $\\hat{\\beta}$", "Standardized 95\\% CI range", 
                          "$\\hat{\\beta}$: ($\\Delta$ men 18-50 x returned to Syria)", "", "  \\emph{p}-value", "  Standarized $\\hat{\\beta}$", "Standardized 95\\% CI range", 
                          "$\\hat{\\beta}$: Family member under siege (yes=1, no=0)", "", "  \\emph{p}-value", "  Standarized $\\hat{\\beta}$", "Standardized 95\\% CI range", 
                          "$\\hat{\\beta}$: (Returned x under siege)", "", "  \\emph{p}-value", "  Standarized $\\hat{\\beta}$", "Standardized 95\\% CI range"
) 


# Create empty vector for control-group means
controlMeans_short1 <- rep("", times = length(RowsNameswSEs_short1))

# Add the control-group means to every 3rd cell
controlMeans_short1[seq(from = 1, to = (numVar_short*5)-4, by = 5) ] <- basicStats_short_control[,"Mean"]


# Define empty results_short1 matrix for tex output
results_short1 <- matrix(nrow = numVar_short*5, ncol = 4)

# add coefficients to tex output
results_short1[seq(from = 1, to = (numVar_short*5)-4, by = 5), ] <- coefs_short

# add CR SEs to text output
results_short1[seq(from = 2, to = (numVar_short*5)-3, by = 5), ] <- paste0("(", sesCR_short, ")")

# add p values  to text output
results_short1[seq(from = 3, to = (numVar_short*5)-2, by = 5), ] <- ifelse(test = p_valuesCR_short_rounded>0.01, yes = p_valuesCR_short_rounded, no = "<0.01")


cohensD_short1_beta <- coefs_short/apply(X = outcomesData_short[, 1:numVar_short], MARGIN = 2, FUN = sd)
results_short1[seq(from = 4, to = (numVar_short*5)-1, by = 5), ] <- abs(round(cohensD_short1_beta, digits = 3))

standardized_CI_range <- (3.92*sesCR_short)/apply(X = outcomesData_short[,1:numVar_short], MARGIN = 2, FUN = sd)
results_short1[seq(from = 5, to = (numVar_short*5), by = 5), ] <- abs(round(standardized_CI_range, digits = 3))



# Bind names, results_short1, and control-group means
results_short1Names <- cbind(RowsNameswSEs_short1, results_short1, controlMeans_short1)


# Add column names

# Add row at the bottom showing whether Controls were included in each model
results_short1Names <- rbind(results_short1Names, c("Covariates", "No", "Yes", "No", "Yes", ""))

colnames(results_short1Names) <- c(" ", "Linear", "Linear w Covariates", "Quadratic", "Quadratic w Covariates", "Control-Group Mean")

#Create duplicate table with unrounded pvalues
results_short1_unrounded <- results_short1
results_short1_unrounded[seq(from = 3, to = (numVar_short*5)-2, by = 5), ] <- p_valuesCR_short



# Output tex table 
print(
  xtable(results_short1Names),
  only.contents = TRUE,
  include.rownames = F,
  include.colnames = F,
  floating = F,
  sanitize.text.function = identity,
  sanitize.colnames.function = identity,
  hline.after = c(
    numVar_short,
    numVar_short * 2,
    numVar_short * 3,
    numVar_short * 4,
    numVar_short * 5
  ),
  file = 'output/results_short1.tex'
)




########## Add appendix table 2

###### For appendix
# Output short2 results table for all models, first five variables
# All models
# columns are models
# rows are point estimate, standard error of the point estimate, p-value, quasi-cohen's d for beta hat, and quasi-cohen's d for beta hat confidence interval range

coefs_short2 <- coefs_long_unrounded[(numVar_short+1):numVar_long, ]
sesCR_short2 <- sesCR_long_unrounded[(numVar_short+1):numVar_long, ]
p_valuesCR_short2 <- p_valuesCR_long[(numVar_short+1):numVar_long, ]

p_valuesCR_short2_rounded <- round(p_valuesCR_short2, digits = 2)

RowsNameswSEs_short2 <- c("$\\hat{\\beta}$: Sum men ages 18 to 50", "", "  \\emph{p}-value", "  Standarized $\\hat{\\beta}$", "Standardized 95\\% CI range", 
                          "$\\hat{\\beta}$: Dangerous work (yes=1, no=0)", "", "  \\emph{p}-value", "  Standarized $\\hat{\\beta}$", "Standardized 95\\% CI range", 
                          "$\\hat{\\beta}$: Went back to syria for money (yes=1, no=0)", "", "  \\emph{p}-value", "  Standarized $\\hat{\\beta}$", "Standardized 95\\% CI range", 
                          "$\\hat{\\beta}$: Money from Syria (yes=1, no=0)", "", "  \\emph{p}-value", "  Standarized $\\hat{\\beta}$", "Standardized 95\\% CI range", 
                          "$\\hat{\\beta}$: $\\Delta$ men among HHs with a return to a war zone", "", "  \\emph{p}-value", "  Standarized $\\hat{\\beta}$", "Standardized 95\\% CI range")



# Create empty vector for control-group means
controlMeans_short2 <- rep("", times = length(RowsNameswSEs_short2))

# Add the control-group means to every 3rd cell
controlMeans_short2[seq(from = 1, to = (numVar_short2*5)-4, by = 5) ] <- basicStats_long_control[(numVar_short+1):numVar_long, "Mean"]


# Define empty results_short2 matrix for tex output
results_short2 <- matrix(nrow = numVar_short2*5, ncol = 4)

# add coefficients to tex output
results_short2[seq(from = 1, to = (numVar_short2*5)-4, by = 5), ] <- coefs_short2

results_short2[seq(from = 1, to = (numVar_short2*5)-4, by = 5), ] <- 
  ifelse(
    round(coefs_short2, digits = 3) == 0,
    yes = formatC(coefs_short2, format = "e", digits = 0),
    no = round(coefs_short2, digits = 3)
  )


# add CR SEs to text output
results_short2[seq(from = 2, to = (numVar_short2*5)-3, by = 5), ] <- paste0("(", round(sesCR_short2, digits = 3), ")")

# add p values  to text output
results_short2[seq(from = 3, to = (numVar_short2*5)-2, by = 5), ] <- ifelse(test = p_valuesCR_short2_rounded>0.01, yes = p_valuesCR_short2_rounded, no = "<0.01")

cohensD_short2_beta <- coefs_short2/apply(X = outcomesData_long[, (numVar_short+1):numVar_long], MARGIN = 2, FUN = sd)
results_short2[seq(from = 4, to = (numVar_short2*5)-1, by = 5), ] <-   abs(round(cohensD_short2_beta, digits = 3))

ifelse(
  abs(round(cohensD_short2_beta, digits = 3)) == 0,
  yes = formatC(abs(cohensD_short2_beta), format = "e", digits = 0),
  no = abs(round(cohensD_short2_beta, digits = 3))
)

standardized_CI_range2 <- (3.92*sesCR_short2)/apply(X = outcomesData_long[, (numVar_short+1):numVar_long], MARGIN = 2, FUN = sd)
results_short2[seq(from = 5, to = (numVar_short2*5), by = 5), ] <- abs(round(standardized_CI_range2, digits = 3))



# Bind names, results_short2, and control-group means
results_short2Names <- cbind(RowsNameswSEs_short2, results_short2, controlMeans_short2)


# Add column names

# Add row at the bottom showing whether Controls were included in each model
results_short2Names <- rbind(results_short2Names, c("Covariates", "No", "Yes", "No", "Yes", ""))

colnames(results_short2Names) <- c(" ", "Linear", "Linear w Covariates", "Quadratic", "Quadratic w Covariates", "Control-Group Mean")


#Create duplicate table with unrounded pvalues
results_short2_unrounded <- results_short2
results_short2_unrounded[seq(from = 3, to = (numVar_short*5)-2, by = 5), ] <- p_valuesCR_short2







# Output tex table 
print(
  xtable(results_short2Names),
  only.contents = TRUE,
  include.rownames = F,
  include.colnames = F,
  floating = F,
  sanitize.text.function = identity,
  sanitize.colnames.function = identity,
  hline.after = c(
    numVar_short2,
    numVar_short2 * 2,
    numVar_short2 * 3,
    numVar_short2 * 4,
    numVar_short2 * 5
  ),
  file = 'output/results_short2.tex'
)


####
# PCA Analysis Results ----------------------------------------------------
#  results_short_PCA.tex 

########## Add appendix table 2

###### For appendix
# Output short_PCA results table for all models, first five variables
# All models
# columns are models
# rows are point estimate, standard error of the point estimate, p-value, quasi-cohen's d for beta hat, and quasi-cohen's d for beta hat confidence interval range


#data[, PCA1  := PCA1]
#data[, PCA2  := PCA2]
#data[, PCA3  := PCA3]
#data[, PCA_other  := PCA_other]

coefs_short_PCA <- coefs_long[(numVar_long+1):(numVar_long+numVar_PCA), ]
sesCR_short_PCA <- sesCR_long[(numVar_long+1):(numVar_long+numVar_PCA), ]
p_valuesCR_short_PCA <- p_valuesCR_long[(numVar_long+1):(numVar_long+numVar_PCA), ]

p_valuesCR_short_PCA_rounded <- round(p_valuesCR_short_PCA, digits = 2)


RowsNameswSEs_short_PCA <- c("$\\hat{\\beta}$: PCA dim. 1", "", "  p-value",     #"  Standarized $\\hat{\\beta}$", "Standardized 95\\% CI range", 
                             "$\\hat{\\beta}$: PCA dim. 2", "", "  p-value",     #"  Standarized $\\hat{\\beta}$", "Standardized 95\\% CI range", 
                             #         "$\\hat{\\beta}$: PCA dim. 3", "", "  p-value",    #"  Standarized $\\hat{\\beta}$", "Standardized 95\\% CI range", 
                             "$\\hat{\\beta}$: PCA dims. 3-8", "", "  p-value"   #, "  Standarized $\\hat{\\beta}$", "Standardized 95\\% CI range"
)


#
# Create empty vector for control-group means

# Add the control-group means to every 3rd cell


# Define empty results_short_PCA matrix for tex output
results_short_PCA <- matrix(nrow = numVar_PCA*3, ncol = 4)

# add coefficients to tex output
results_short_PCA[seq(from = 1, to = (numVar_PCA*3)-2, by = 3), ] <- coefs_short_PCA

# add CR SEs to text output
results_short_PCA[seq(from = 2, to = (numVar_PCA*3)-1, by = 3), ] <- paste0("(", sesCR_short_PCA, ")")

# add p values  to text output
results_short_PCA[seq(from = 3, to = (numVar_PCA*3), by = 3), ] <- ifelse(test = p_valuesCR_short_PCA_rounded>0.01, yes = p_valuesCR_short_PCA_rounded, no = "<0.01")


# Load PCA values
eig1 <- as.numeric(read.csv("output/PCA_eig1.csv"))
eig2 <- as.numeric(read.csv("output/PCA_eig2.csv"))
eig_others <- as.numeric(read.csv("output/PCA_eig_others.csv"))
eigs <- c(eig1, eig2, eig_others)

eigs_wSpaces <- rep("", times = numVar_PCA*3)
eigs_wSpaces[seq(from = 1, to = (numVar_PCA*3)-2, by = 3) ] <- round(eigs, digits = 2)


# Bind names, results_short_PCA
results_short_PCANames <- cbind(RowsNameswSEs_short_PCA, results_short_PCA, eigs_wSpaces)


# Add column names

# Add row at the bottom showing whether Controls were included in each model
results_short_PCANames <- rbind(results_short_PCANames, c("Controls", "No", "Yes", "No", "Yes", ""))

#
colnames(results_short_PCANames) <- c(" ", "Linear", "Linear w Controls", "Quadratic", "Quadratic w Controls", "\\% Var. Exp.")

# Output tex table 
print(
  xtable(results_short_PCANames),
  only.contents = TRUE,
  include.rownames = F,
  include.colnames = F,
  floating = F,
  sanitize.text.function = identity,
  sanitize.colnames.function = identity,
  hline.after = c(3, 6, 9),
  file = 'output/results_short_PCA.tex'
)













####
# Tables for Main Paper ---------------------------------------------------

#   linear_withControls_beta_SE_meanControl_pvalue.tex 
## Now produce short table with just linear w covariates

RowsNameswSEs_short_linearC <- paste("(", sprintf("%.3f", sesCR_short[, "Linear_C"]), ")", sep="") 

cohensD_short_beta_linear_C <- coefs_short[, "Linear_C"]/apply(X = outcomesData_short[Treatment==0, 1:numVar_short], MARGIN = 2, FUN = sd)

standardized_SE_linearWcovs_holder <- sesCR_short[, "Linear_C"]/apply(X = outcomesData_short[Treatment==0, 1:numVar_short], MARGIN = 2, FUN = sd)
standardized_SE_linearWcovs_holder2 <- round(standardized_SE_linearWcovs_holder, digits = 3)
standardized_SE_linearWcovs <- paste0("[", standardized_SE_linearWcovs_holder2, "]") 


linear_c_table <- rbind(coefs_short[, "Linear_C"], 
                        RowsNameswSEs_short_linearC,
                        basicStats_short_control[,"Mean"], 
                        round(p_valuesCR_short[, "Linear_C"], digits=3),
                        round(cohensD_short_beta_linear_C, digits = 3),
                        standardized_SE_linearWcovs
)

# high estimate of treatment group outcome
as.numeric(linear_c_table[3,])+(as.numeric(linear_c_table[1,])+1.96*sesCR_short[, "Linear_C"])
# low estimate of treatment group outcome
as.numeric(linear_c_table[3,])+(as.numeric(linear_c_table[1,])-1.96*sesCR_short[, "Linear_C"])

rownames(linear_c_table) <- c("$\\hat{\\beta}$",
                              "",
                              "Control-group mean",
                              "\\emph{p}-value",
                              "Standardized $\\hat{\\beta}$",
                              " " #"Standardized 95\\% CI range"
)
colnames(linear_c_table) <- c("1. Change in number of men ages 18-50",
                              "2. A household member returned to Syria", 
                              "3. Change in number of men 18-50 when someone returned to Syria",
                              "4. Family member in an active war zone", 
                              "5. Household member returned to an active war zone")

# Output tex table 
print(xtable(linear_c_table), only.contents=TRUE, include.rownames=T, 
      include.colnames=T, floating=F, sanitize.rownames.function = identity,
      file = 'output/linear_withControls_beta_SE_meanControl_pvalue.tex')


##################################
# Plots for Main Paper ----------------------------------------------------

## Output RD plots
pdf("output/rdplot_linear_main_outcomes.pdf", width=6, height=4) 
par(mfrow=c(2,3))

rdplot(y = outcomesData_short[, 1], x = outcomesData_short $WINTELEV_oct13_rescaled, c = 0, p = 1, binselect = "esmvpr", x.label = "Var. 1", y.label = "")
rdplot(y = outcomesData_short[, 2], x = outcomesData_short $WINTELEV_oct13_rescaled, c = 0, p = 1, binselect = "esmvpr", x.label = "Var. 2", y.label = "")
rdplot(y = outcomesData_short[, 3], x = outcomesData_short $WINTELEV_oct13_rescaled, c = 0, p = 1, binselect = "esmvpr", x.label = "Var. 3", y.label = "")
rdplot(y = outcomesData_short[, 4], x = outcomesData_short $WINTELEV_oct13_rescaled, c = 0, p = 1, binselect = "esmvpr", x.label = "Var. 4", y.label = "")
rdplot(y = outcomesData_short[, 5], x = outcomesData_short $WINTELEV_oct13_rescaled, c = 0, p = 1, binselect = "esmvpr", x.label = "Var. 5", y.label = "")

dev.off()



require(dplyr)

var1_bin_means <-    as.data.frame(outcomesData_long %>%
                                     group_by(WINTELEV_oct13_rescaled) %>%
                                     summarise(mean_by_alt = mean(DiD18_50), n = n(), n_scale_holder = log(n()) ))

var1_bin_means$n_scale <-   var1_bin_means$n_scale_holder/max(var1_bin_means$n_scale_holder)


var2_bin_means <-  as.data.frame(outcomesData_long %>%
                                   group_by(WINTELEV_oct13_rescaled) %>%
                                   summarise(mean_by_alt = mean(`Q64>0`), n = n(), n_scale_holder = log(n()) ))

var2_bin_means$n_scale <-   var2_bin_means$n_scale_holder/max(var2_bin_means$n_scale_holder)


var3_bin_means <-  as.data.frame(outcomesData_long %>%
                                   group_by(WINTELEV_oct13_rescaled) %>%
                                   summarise(mean_by_alt = mean(`I(rescaled_men_18_50xQ64_0_1)`), n = n(), n_scale_holder = log(n()) ))

var3_bin_means$n_scale <-   var3_bin_means$n_scale_holder/max(var3_bin_means$n_scale_holder)


var4_bin_means <-  as.data.frame(outcomesData_long %>%
                                   group_by(WINTELEV_oct13_rescaled) %>%
                                   summarise(mean_by_alt = mean(`66b1`), n = n(), n_scale_holder = log(n()) ))

var4_bin_means$n_scale <-   var4_bin_means$n_scale_holder/max(var4_bin_means$n_scale_holder)


var5_bin_means <-  as.data.frame(outcomesData_long %>%
                                   group_by(WINTELEV_oct13_rescaled) %>%
                                   summarise(mean_by_alt = mean(returnedXundersiege), n = n(), n_scale_holder = log(n()) ))

var5_bin_means$n_scale <-   var5_bin_means$n_scale_holder/max(var5_bin_means$n_scale_holder)




##############################################
# Output plot for models with covariate adjustment -------------------------------
pdf("output/rdplot_linear_main_outcomes_w_covs.pdf", width=9, height=6) 

par(mfrow=c(2,3), oma=c(0,0,2,0))

#c <- 0
x <- outcomesData_long[, "WINTELEV_oct13_rescaled"]
#y0 <- round(x*2.5 + rnorm(100))
Z <- outcomesData_long[, "Treatment"]

# var 1
y <- outcomesData_long[, 1]

model <- lm(y ~ Z*x   +    outcomesData_long$age0to4 +
              outcomesData_long$age5to12 +
              outcomesData_long$age13to17 +
              outcomesData_long$age18to59 +
              outcomesData_long$age60plus +
              outcomesData_long$Oct2013Pop +
              outcomesData_long$disabled +
              outcomesData_long$origin_gov_0 +
              outcomesData_long$origin_gov_4 +
              outcomesData_long$origin_gov_6 +
              outcomesData_long$origin_gov_10 +
              outcomesData_long$origin_gov_11 +
              outcomesData_long$educ_le_priinc +
              outcomesData_long$educ_pricom +
              outcomesData_long$educ_middle +
              outcomesData_long$educ_sec_tech_uni +
              outcomesData_long$age_head +
              outcomesData_long$monthinleb +
              outcomesData_long$base_relatives +
              outcomesData_long$base_friends )

a <- model$coefficients['(Intercept)']         
betaZ <- model$coefficients['Z']           
betaX <- model$coefficients['x']         
betaZX <- model$coefficients['Z:x'] 


# means   plot(x = var1_bin_means$WINTELEV_oct13_rescaled, y = var1_bin_means$mean_by_alt, xlab = "Altitude", ylab = "Var. 1 mean by altitude",col = alpha('black', 0.2), pch=16, xlim = c(-50,50))
plot(x = x, y = y, xlab = "Altitude", ylab = "Var. 1 by household",col = alpha('black', var1_bin_means$n_scale), pch=16, xlim = c(-50,50))

lines(x = c(min(x), -0.000001), y = c(a + min(x)*betaX, a + -0.000001*betaX), col="red") # regression line (y~x) 
lines(x = c(0, max(x)), y = c(a + betaZ, a + betaZ + max(x)*betaX + max(x)*betaZX), col="red") # regression line (y~x) 



# var 2
y <- outcomesData_long[, 2]

model <- lm(y ~ Z*x   +    outcomesData_long$age0to4 +
              outcomesData_long$age5to12 +
              outcomesData_long$age13to17 +
              outcomesData_long$age18to59 +
              outcomesData_long$age60plus +
              outcomesData_long$Oct2013Pop +
              outcomesData_long$disabled +
              outcomesData_long$origin_gov_0 +
              outcomesData_long$origin_gov_4 +
              outcomesData_long$origin_gov_6 +
              outcomesData_long$origin_gov_10 +
              outcomesData_long$origin_gov_11 +
              outcomesData_long$educ_le_priinc +
              outcomesData_long$educ_pricom +
              outcomesData_long$educ_middle +
              outcomesData_long$educ_sec_tech_uni +
              outcomesData_long$age_head +
              outcomesData_long$monthinleb +
              outcomesData_long$base_relatives +
              outcomesData_long$base_friends )

a <- model$coefficients['(Intercept)']         
betaZ <- model$coefficients['Z']           
betaX <- model$coefficients['x']         
betaZX <- model$coefficients['Z:x'] 

plot(x = var2_bin_means$WINTELEV_oct13_rescaled, y = var2_bin_means$mean_by_alt, xlab = "Altitude", ylab = "Var. 2 mean by altitude", col = alpha('black', var2_bin_means$n_scale), pch=16, xlim = c(-50,50))

lines(x = c(min(x), -0.000001), y = c(a + min(x)*betaX, a + -0.000001*betaX), col="red") # regression line (y~x) 
lines(x = c(0, max(x)), y = c(a + betaZ, a + betaZ + max(x)*betaX + max(x)*betaZX), col="red") # regression line (y~x) 




# var 3
y <- outcomesData_long[, 3]

model <- lm(y ~ Z*x   +    outcomesData_long$age0to4 +
              outcomesData_long$age5to12 +
              outcomesData_long$age13to17 +
              outcomesData_long$age18to59 +
              outcomesData_long$age60plus +
              outcomesData_long$Oct2013Pop +
              outcomesData_long$disabled +
              outcomesData_long$origin_gov_0 +
              outcomesData_long$origin_gov_4 +
              outcomesData_long$origin_gov_6 +
              outcomesData_long$origin_gov_10 +
              outcomesData_long$origin_gov_11 +
              outcomesData_long$educ_le_priinc +
              outcomesData_long$educ_pricom +
              outcomesData_long$educ_middle +
              outcomesData_long$educ_sec_tech_uni +
              outcomesData_long$age_head +
              outcomesData_long$monthinleb +
              outcomesData_long$base_relatives +
              outcomesData_long$base_friends )

a <- model$coefficients['(Intercept)']         
betaZ <- model$coefficients['Z']           
betaX <- model$coefficients['x']         
betaZX <- model$coefficients['Z:x'] 

# means plot(x = var3_bin_means$WINTELEV_oct13_rescaled, y = var3_bin_means$mean_by_alt, xlab = "Altitude", ylab = "Var. 3 mean by altitude",col = alpha('black', 0.2), pch=16, xlim = c(-50,50), ylim = c(-.02, .12))
plot(x = x, y = y, xlab = "Altitude", ylab = "Var. 3 by household",col = alpha('black', var3_bin_means$n_scale), pch=16, xlim = c(-50,50) )   #, ylim = c(-.02, .12))


lines(x = c(min(x), -0.000001), y = c(a + min(x)*betaX, (as.numeric(a) + -0.000001*betaX)), col="red") # regression line (y~x) 
lines(x = c(0, max(x)), y = c(a + betaZ, a + betaZ + max(x)*betaX + max(x)*betaZX), col="red") # regression line (y~x) 


lines(x = c(-50, -0.000001), y = c(0.08260675, 0.10827), col="red")

# var 4
y <- outcomesData_long[, 4]

model <- lm(y ~ Z*x   +    outcomesData_long$age0to4 +
              outcomesData_long$age5to12 +
              outcomesData_long$age13to17 +
              outcomesData_long$age18to59 +
              outcomesData_long$age60plus +
              outcomesData_long$Oct2013Pop +
              outcomesData_long$disabled +
              outcomesData_long$origin_gov_0 +
              outcomesData_long$origin_gov_4 +
              outcomesData_long$origin_gov_6 +
              outcomesData_long$origin_gov_10 +
              outcomesData_long$origin_gov_11 +
              outcomesData_long$educ_le_priinc +
              outcomesData_long$educ_pricom +
              outcomesData_long$educ_middle +
              outcomesData_long$educ_sec_tech_uni +
              outcomesData_long$age_head +
              outcomesData_long$monthinleb +
              outcomesData_long$base_relatives +
              outcomesData_long$base_friends )

a <- model$coefficients['(Intercept)']         
betaZ <- model$coefficients['Z']           
betaX <- model$coefficients['x']         
betaZX <- model$coefficients['Z:x'] 

plot(x = var4_bin_means$WINTELEV_oct13_rescaled, y = var4_bin_means$mean_by_alt, xlab = "Altitude", ylab = "Var. 4 mean by altitude",col = alpha('black', var4_bin_means$n_scale), pch=16, xlim = c(-50,50))

lines(x = c(min(x), -0.000001), y = c(a + min(x)*betaX, a + -0.000001*betaX), col="red") # regression line (y~x) 
lines(x = c(0, max(x)), y = c(a + betaZ, a + betaZ + max(x)*betaX + max(x)*betaZX), col="red") # regression line (y~x) 



# var 5
y <- outcomesData_long[, 5]

model <- lm(y ~ Z*x   +    outcomesData_long$age0to4 +
              outcomesData_long$age5to12 +
              outcomesData_long$age13to17 +
              outcomesData_long$age18to59 +
              outcomesData_long$age60plus +
              outcomesData_long$Oct2013Pop +
              outcomesData_long$disabled +
              outcomesData_long$origin_gov_0 +
              outcomesData_long$origin_gov_4 +
              outcomesData_long$origin_gov_6 +
              outcomesData_long$origin_gov_10 +
              outcomesData_long$origin_gov_11 +
              outcomesData_long$educ_le_priinc +
              outcomesData_long$educ_pricom +
              outcomesData_long$educ_middle +
              outcomesData_long$educ_sec_tech_uni +
              outcomesData_long$age_head +
              outcomesData_long$monthinleb +
              outcomesData_long$base_relatives +
              outcomesData_long$base_friends )

a <- model$coefficients['(Intercept)']         
betaZ <- model$coefficients['Z']           
betaX <- model$coefficients['x']         
betaZX <- model$coefficients['Z:x'] 

plot(x = var5_bin_means$WINTELEV_oct13_rescaled, y = var5_bin_means$mean_by_alt, xlab = "Altitude", ylab = "Var. 5 mean by altitude",col = alpha('black', var5_bin_means$n_scale), pch=16, xlim = c(-50,50))

lines(x = c(min(x), -0.000001), y = c(a + min(x)*betaX, a + -0.000001*betaX), col="red") # regression line (y~x) 
lines(x = c(0, max(x)), y = c(a + betaZ, a + betaZ + max(x)*betaX + max(x)*betaZX), col="red") # regression line (y~x) 


frame()
text(x = 0.28, y = 0.95, "Var. 1: Change in number \n           of men ages 18-50")  
text(x = 0.32, y = 0.725, "Var. 2: A household member  \n          returned to Syria (0/1)")
text(x = 0.46, y = 0.5, "Var. 3: Change in number of men 18-50 \n         when someone returned to Syria")
text(x = 0.3, y = 0.275, "Var. 4: A family member in  \n               an active war zone (0/1)")
text(x = 0.3, y = 0.052, "Var. 5: Household member \n                                   returned to an active war zone (0/1)")


dev.off()



# Output Numbers for Results Interpretation Section


####### SIGN
# Looking at the sign of our point estimates, 
#  23 of 40 treatment effect estimates have a sign pointing toward an decrease in mobilization. 
#    DECREASE SIGN
#   CHECK SIGNS and MEANINGS
write( sum(results_short1[1, ]>0) + sum(results_short1[6, ]<0) + sum(results_short1[11, ]>0) + sum(results_short1[16, ]<0) + sum(results_short1[21, ]<0)+
  sum(results_short2[1, ]>0) + sum(results_short2[6, ]<0) + sum(results_short2[11, ]<0) + sum(results_short2[16, ]<0) + sum(results_short2[21, ]>0),
file = "output/estimates_sign_decrease.txt")


#  Among these 23, 
#   XX 8 estimates have \emph{p} values less than 0.05, 
# sum(results_short1[1, ]<0 & results_short1[3, ]<0.05)
write( sum(results_short1[1, ]>0 & results_short1[3, ]<0.05) + sum(results_short1[6, ]<0 & results_short1[8, ]<0.05) + 
  sum(results_short1[11, ]>0 & results_short1[13, ]<0.05) + sum(results_short1[16, ]<0 & results_short1[18, ]<0.05) + 
  sum(results_short1[21, ]<0 & results_short1[23, ]<0.05) +
  sum(results_short2[1, ]>0 & results_short2[3, ]<0.05) + sum(results_short2[6, ]<0 & results_short2[8, ]<0.05) + 
  sum(results_short2[11, ]<0 & results_short2[13, ]<0.05) + sum(results_short2[16, ]<0 & results_short2[18, ]<0.05) + 
  sum(results_short2[21, ]>0 & results_short2[23, ]<0.05),
file = 'output/estimates_sign_decrease_significant.txt')

#  10 have  standardized coefficient sizes greater than 0.2. 
# sum(results_short1[1, ]<0 & results_short1[4, ]>0.2)
write( sum(results_short1[1, ]>0 & results_short1[4, ]>0.2) + sum(results_short1[6, ]<0 & results_short1[9, ]>0.2) + 
  sum(results_short1[11, ]>0 & results_short1[14, ]>0.2) + sum(results_short1[16, ]<0 & results_short1[19, ]>0.2) + 
  sum(results_short1[21, ]<0 & results_short1[24, ]>0.2)+
  sum(results_short2[1, ]>0 & results_short2[4, ]>0.2) + sum(results_short2[6, ]<0 & results_short2[9, ]>0.2) + 
  sum(results_short2[11, ]<0 & results_short2[14, ]>0.2) + sum(results_short2[16, ]<0 & results_short2[19, ]>0.2) + 
  sum(results_short2[21, ]>0 & results_short2[24, ]>0.2),
file = 'output/estimates_sign_decrease_not_small.txt')

#   6 treatment effect estimates both have \emph{p} values less than 0.05 and  standardized coefficient sizes greater than 0.2. 
#  sum(results_short1[1, ]<0 & results_short1[3, ]<0.05 & results_short1[4, ]>0.2)
write(sum(results_short1[1, ]>0 & results_short1[3, ]<0.05 & results_short1[4, ]>0.2) + 
  sum(results_short1[6, ]<0 & results_short1[8, ]<0.05 & results_short1[9, ]>0.2) + 
  sum(results_short1[11, ]>0 & results_short1[13, ]<0.05 & results_short1[14, ]>0.2) + 
  sum(results_short1[16, ]<0 & results_short1[18, ]<0.05 & results_short1[19, ]>0.2) + 
  sum(results_short1[21, ]<0 & results_short1[23, ]<0.05 & results_short1[24, ]>0.2)+
  sum(results_short2[1, ]>0 & results_short2[3, ]<0.05 & results_short2[4, ]>0.2) + 
  sum(results_short2[6, ]<0 & results_short2[8, ]<0.05 & results_short2[9, ]>0.2) + 
  sum(results_short2[11, ]<0 & results_short2[13, ]<0.05 & results_short2[14, ]>0.2) + 
  sum(results_short2[16, ]<0 & results_short2[18, ]<0.05 & results_short2[19, ]>0.2) + 
  sum(results_short2[21, ]>0 & results_short2[23, ]<0.05 & results_short2[24, ]>0.2),
file = 'output/estimates_sign_decrease_significant_not_small.txt')


#  These 6 estimates are entirely for two outcomes variables: 
# whether a household member undertook dangerous work in the past six months, and 
# the interaction of whether a household member returned and whether a family member is living in a siege zone. 

###  SIGN pointing toward increase
#  Among the  17  of 40 treatment effect estimates that have a sign pointing toward an increase in mobilization, 
# almost all are well identified zeros.  
#    INCREASE SIGN
#   CHECK SIGNS and MEANINGS
write( sum(results_short1[1, ]<0) + sum(results_short1[6, ]>0) + sum(results_short1[11, ]<0) + sum(results_short1[16, ]>0) + sum(results_short1[21, ]>0) +
  sum(results_short2[1, ]<0) + sum(results_short2[6, ]>0) + sum(results_short2[11, ]>0) + sum(results_short2[16, ]>0) + sum(results_short2[21, ]<0),
  file = 'output/estimates_sign_increase.txt')



#  Among these 17, 
#only 
#  1 estimate has a \emph{p} value less than 0.05, 
#  sum(results_short1[1, ]>0 & results_short1[3, ]<0.05)
write( sum(results_short1[1, ]<0 & results_short1[3, ]<0.05) + sum(results_short1[6, ]>0 & results_short1[8, ]<0.05) + 
  sum(results_short1[11, ]<0 & results_short1[13, ]<0.05) + sum(results_short1[16, ]>0 & results_short1[18, ]<0.05) + 
  sum(results_short1[21, ]>0 & results_short1[23, ]<0.05) +
  sum(results_short2[1, ]<0 & results_short2[3, ]<0.05) + sum(results_short2[6, ]>0 & results_short2[8, ]<0.05) + 
  sum(results_short2[11, ]>0 & results_short2[13, ]<0.05) + sum(results_short2[16, ]>0 & results_short2[18, ]<0.05) + 
  sum(results_short2[21, ]<0 & results_short2[23, ]<0.05),
file = 'output/estimates_sign_increase_significant.txt')

# % of treatment effect estimates where to sign points toward an increase in mobilization,
write(( sum(results_short1[1, ]<0) + sum(results_short1[6, ]>0) + sum(results_short1[11, ]<0) + 
          sum(results_short1[16, ]>0) + sum(results_short1[21, ]>0) +
          sum(results_short2[1, ]<0) + sum(results_short2[6, ]>0) + sum(results_short2[11, ]>0) +
          sum(results_short2[16, ]>0) + sum(results_short2[21, ]<0))/40*100,
      file = "output/share_increase.txt")


#  and all 17 have standardized coefficient sizes less than 0.2.
#    sum(results_short1[1, ]>0 & results_short1[4, ]>0.2)
sum(results_short1[1, ]<0 & results_short1[4, ]>0.2) + sum(results_short1[6, ]>0 & results_short1[9, ]>0.2) + 
  sum(results_short1[11, ]<0 & results_short1[14, ]>0.2) + sum(results_short1[16, ]>0 & results_short1[19, ]>0.2) + 
  sum(results_short1[21, ]>0 & results_short1[24, ]>0.2) +
  sum(results_short2[1, ]<0 & results_short2[4, ]>0.2) + sum(results_short2[6, ]>0 & results_short2[9, ]>0.2) + 
  sum(results_short2[11, ]>0 & results_short2[14, ]>0.2) + sum(results_short2[16, ]>0 & results_short2[19, ]>0.2) + 
  sum(results_short2[21, ]<0 & results_short2[24, ]>0.2)



#  That is, among the 42.5\% of treatment effect estimates where to sign points toward an increase in mobilization, 
#  they are all substantively small and only one is statistically distinguishable from zero.

###### p Values
# Looking at \emph{p} values of our point estimates, 
#  8 of 40 treatment effect estimates have \emph{p} values less than 0.05, 
write( sum(results_short1[3, ]<0.05) + sum(results_short1[8, ]<0.05) + sum(results_short1[13, ]<0.05) + sum(results_short1[18, ]<0.05) + sum(results_short1[23, ]<0.05) +
  sum(results_short2[3, ]<0.05) + sum(results_short2[8, ]<0.05) + sum(results_short2[13, ]<0.05) + sum(results_short2[18, ]<0.05) + sum(results_short2[23, ]<0.05),
file = 'output/estimates_significant.txt')


#  Among the 20\% treatment effect estimates that are statistically distinguishable from zero, 
#    all but two point toward a decrease in mobilization.
write( (sum(results_short1[3, ]<0.05) + sum(results_short1[8, ]<0.05) + sum(results_short1[13, ]<0.05) + sum(results_short1[18, ]<0.05) + sum(results_short1[23, ]<0.05) +
         sum(results_short2[3, ]<0.05) + sum(results_short2[8, ]<0.05) + sum(results_short2[13, ]<0.05) + sum(results_short2[18, ]<0.05) + sum(results_short2[23, ]<0.05))/40*100,
      file = "output/share_significant.txt")

#   6 of these 8 estimates point toward a negative treatment effect on mobilization. 
#   sum(results_short1[3, ]>0 & results_short1[1, ]<0)
write( sum(results_short1[3, ]<0.05 & results_short1[1, ]>0) + sum(results_short1[8, ]<0.05 & results_short1[6, ]<0) + 
  sum(results_short1[13, ]<0.05 & results_short1[11, ]>0) + sum(results_short1[18, ]<0.05 & results_short1[16, ]<0) + 
  sum(results_short1[23, ]<0.05 & results_short1[21, ]<0) +
  sum(results_short2[3, ]<0.05 & results_short2[1, ]>0) + sum(results_short2[8, ]<0.05 & results_short2[6, ]<0) + 
  sum(results_short2[13, ]<0.05 & results_short2[11, ]<0) + sum(results_short2[18, ]<0.05 & results_short2[16, ]<0) + 
  sum(results_short2[23, ]<0.05 & results_short2[21, ]>0),
  file = 'output/estimates_significant_sign_decrease.txt')






##### Standardized coefficient size, 
# 10 of 40 treatment effect estimates have a standardized treatment effect estimate larger than 0.2 
write( sum(results_short1[4, ]>0.2) + sum(results_short1[9, ]>0.2) + sum(results_short1[14, ]>0.2) + sum(results_short1[19, ]>0.2) + sum(results_short1[24, ]>0.2)+ 
  sum(results_short2[4, ]>0.2) + sum(results_short2[9, ]>0.2) + sum(results_short2[14, ]>0.2) + sum(results_short2[19, ]>0.2) + sum(results_short2[24, ]>0.2),
file = 'output/estimates_not_small.txt')


# among these 10 regressions,  all  estimates point in the direction of a decrease in mobilization. 
#  sum(results_short1[4, ]>0.2 & results_short1[1, ]<0)
write( sum(results_short1[4, ]>0.2 & results_short1[1, ]>0) + sum(results_short1[9, ]>0.2 & results_short1[6, ]<0) + 
  sum(results_short1[14, ]>0.2 & results_short1[11, ]>0) + sum(results_short1[19, ]>0.2 & results_short1[16, ]<0) + 
  sum(results_short1[24, ]>0.2 & results_short1[21, ]<0)+ 
  sum(results_short2[4, ]>0.2 & results_short2[1, ]>0) + sum(results_short2[9, ]>0.2 & results_short2[6, ]<0) + 
  sum(results_short2[14, ]>0.2 & results_short2[11, ]<0) + sum(results_short2[19, ]>0.2 & results_short2[16, ]<0) +
  sum(results_short2[24, ]>0.2 & results_short2[21, ]>0),
file = 'output/estimates_not_small_sign_decrease.txt')

# Among the 25\% of treatment effect estimates that are not substantively small, 
write( (sum(results_short1[4, ]>0.2) + sum(results_short1[9, ]>0.2) + sum(results_short1[14, ]>0.2) + 
         sum(results_short1[19, ]>0.2) + sum(results_short1[24, ]>0.2)+ 
         sum(results_short2[4, ]>0.2) + sum(results_short2[9, ]>0.2) + sum(results_short2[14, ]>0.2) + 
         sum(results_short2[19, ]>0.2) + sum(results_short2[24, ]>0.2))/40*100,
       file = 'output/share_not_small.txt')




